\label{sec:app}

As mentioned in Section~\ref{sec:model}, the model is built upon KCL,
therefore, the current equations for each cross-point can be set following
\begin{equation}\label{equ:KCL0}
 {\Sigma}_{I=1}^kI_k=0.
\end{equation}
For the sake of brevity, we assume that the wordline voltage drivers are
only located at the edge of $V_{W1} \sim V_{Wm}$ and bitline multiplexors
or are located at the edge of $V_{B1} \sim V_{Bn}$. Points located at the
other two edges are left floating.

First of all, for normal points which are located inside the memory array,
the KCL equations take the form of
\begin{equation}\label{equ:KCL1}
R_l^{-1}V_{i,j-1} -(2R_l^{-1}+R_{i,j}^{-1})V_{i,j}+ R_l^{-1}V_{i,j+1}+R_{i,j}^{-1}V'_{i,j}=0,
\end{equation}
for the node at wordline layer and
\begin{equation}\label{equ:KCL2}
R_l^{-1}V'_{i-1,j} -(2R_l^{-1}+R_{i,j}^{-1})V'_{i,j}+ R_l^{-1}V'_{i+1,j}+R_{i,j}^{-1}V_{i,j}=0,
\end{equation}
for all of the nodes with $1<i<m$ and $1<j<n$ in a $m \times n$ array.

For all of the points $V_{i,1}$ ($1\leq i\leq m$) according to different
write schemes, they can be connected to the voltage driver $V_{Wi}$ (as
activated points) or left floating (as floating points). For activated
points, we have
\begin{equation}\label{equ:KCL3}
 -(R_v^{-1}+R_l^{-1}+R_{i,1}^{-1})V_{i,1}+ R_l^{-1}V_{i,2}+R_{i,1}^{-1}V'_{i,1}=-R_v^{-1}V_{Wi},
\end{equation}
and for floating points, we have
\begin{equation}\label{equ:KCL4}
 -(R_l^{-1}+R_{i,1}^{-1})V_{i,1}+ R_l^{-1}V_{i,2}+R_{i,j}^{-1}V'_{i,1}=0.
\end{equation}

Similarly, for the points of $V'_{1,j}$ ($1\leq j\leq n$), the KCL
equations take the form of
\begin{equation}\label{equ:KCL5}
 -(R_s^{-1}+R_l^{-1}+R_{1,j}^{-1})V'_{1,j}+ R_l^{-1}V'_{2,j}+R_{1,j}^{-1}V_{1,j}=-R_s^{-1}V_{Bj},
\end{equation}
for activated points and
\begin{equation}\label{equ:KCL6}
 -(R_l^{-1}+R_{1,j}^{-1})V'_{1,j}+ R_l^{-1}V'_{2,j}+R_{1,j}^{-1}V_{1,j}=0.
\end{equation}
for floating points.

Finally, all of the other points at $V_{i,n}$ and $V'_{m,j}$ ($1\leq i\leq
m, 1\leq j\leq n$) are floating points and have the form of
\begin{equation}\label{equ:KCL7}
 -(R_l^{-1}+R_{i,n}^{-1})V_{i,n}+ R_l^{-1}V_{i,n-1}+R_{i,n}^{-1}V'_{i,n}=0,
\end{equation}
\begin{equation}\label{equ:KCL8}
 -(R_l^{-1}+R_{m,j}^{-1})V'_{m,j}+ R_l^{-1}V'_{m-1,j}+R_{m,j}^{-1}V_{m,j}=0.
\end{equation}

Then, for clarity, a ${2mn\times 1}$ vector ${V}$ is defined to represent
all of the variables in the KCL equations:
\begin{equation}\label{equ:V1}
{V}=[{V_1}^T,{V_2}^T...{V_m}^T,{V'_1}^T,{V'_2}^T...{V'_m}^T]^T,
\end{equation}
where,
%\begin{equation}\label{equ:V2}
%{V_i} = [V_{i,1},V_{i,2}...V_{i,n}]^T,\\
%\end{equation}
%\begin{equation}\label{equ:V3}
%{V'_i} = [V'_{i,1},V'_{i,2}...V'_{i,n}]^T,
%\end{equation}
\begin{equation}\label{equ:V2}
{V_i} = [V_{i,1},V_{i,2}...V_{i,n}]^T,~~{V'_i} = [V'_{i,1},V'_{i,2}...V'_{i,n}]^T,
\end{equation}
for $i=1,2...m$. Then all of the KCL equations can be considered as a
system of linear equations, which has the form
\begin{equation}\label{equ:matrix}
A\cdot V = C.
\end{equation}
$A$ is a ${2mn\times{2mn}}$ coefficient matrix, which is determined by
Equations(\ref{equ:KCL1})-(\ref{equ:KCL8}). $C$ is a ${2mn\times{1}}$
vector, containing the constant terms of these equations. Obviously, the
KCL equation for each point has a relatively simple structure and they are
similar to each other. Thus, the linear equation system has a fixed format
and simple structure, which is easy to establish and adjust according to
different design schemes and cell parameters. Moreover, matrix $A$ is
populated primarily with zeros and can be saved as a sparse matrix, which
will further reduce the storage cost during the computation.

The characteristics of the linear system can be summarized as
\begin{enumerate}
  \item As shown in Equation~(\ref{equ:blockedmatrix}), the
      coefficient matrix $A$ can be further partitioned into four
      subblocks
    \begin{equation}\label{equ:blockedmatrix}
        \mathbf{A} = \left[
        \begin{array}{cc}
            A1 & A2  \\
            A3 & A4  \\
        \end{array} \right].
    \end{equation}
All of these subblocks have the same size of $mn\times mn$. Subblock
$A2$ and $A3$ are diagonal matrixes and have the same form of
%\begin{equation}\label{equ:A2A3}
%        A2 = A3 = \left[
%        \begin{array}{cccc}
%            R_{1,1}^{-1}    & 0             & \ldots    & 0 \\
%            0               & R_{2,2}^{-1}  & \ddots    & \vdots  \\
%            \vdots          & \ddots        & \ddots    & 0 \\
%            0               & \ldots        & 0         & R_{mn,mn}^{-1} \\
%        \end{array} \right].
%    \end{equation}
\begin{equation}\label{equ:A2A3}
        A2 = A3 = \left[
        \begin{array}{ccccccc}
            R_{1,1}^{-1}    & 0             & 0             & ~         & \ldots    & ~             & 0    \\
            0               & \ddots        & 0             & ~         & ~         & ~             & ~    \\
            0               & 0             & R_{1,n}^{-1}  & ~         & \ddots    & ~             & \vdots    \\
            ~               & ~             & ~             & \ddots    & ~         & ~             & ~    \\
            \vdots          & ~        & \ddots             & ~         & R_{m,1}^{-1}    & 0       & 0   \\
            ~               & ~             & ~             & ~         & 0         & \ddots        & 0    \\
            0               & 0             & \ldots              & ~         & 0    & 0             & R_{m,n}^{-1}    \\

        \end{array} \right].
\end{equation}
$A2$ and $A3$ do not change their values with different schemes.

However, $A1$ and $A4$ are a little more complex than $A2$ and $A3$.
$A1$ is a tridiagonal matrix and has nonzero elements only located in
the main diagonal, and the first line below and above the diagonal.
Similarly, $A4$ is a special tridiagonal matrix, which has nonzero
elements in the main diagonal, and the $n^{th}$ line below and above
the diagonal, where $n$ is the number of bitlines in the cross-point
model. The coefficients of normal points in $A1$ and $A4$ can be
easily derived from Equation (\ref{equ:KCL1}) and (\ref{equ:KCL2}) as

\begin{equation}\label{equ:A12}
\left\{
    \begin{array}{ll}
    A1(k,k-1) = R_l^{-1}\\
    A1(k,k)~~~~~ = -(2R_l^{-1}+R_{i,j}^{-1})~~~~,\\
    A1(k,k+1) = R_l^{-1}\\
    \end{array} \right.
    \end{equation}

and
\begin{equation}\label{equ:A122}
\left\{
    \begin{array}{ll}
    A4(k,k-n) = R_l^{-1}\\
    A4(k,k)~~~~~ = -(2R_l^{-1}+R_{i,j}^{-1})~~~~,\\
    A4(k,k+n) = R_l^{-1}\\
    \end{array} \right.
    \end{equation}
where $ k=(i-1)n+j$, $1<i<m$, and $1<j<n$.

However, since the edge conditions vary with different program
schemes, the coefficients of points at the edge of the array may be
different. Clearly, the four edges shown in Figure~\ref{fig:modeling}
correspond to different coefficients in $A1$ and $A4$. Therefore,
after initializing matrix $A1$ and $A4$ following
Equations(\ref{equ:A12}) and (\ref{equ:A122}), the coefficients
related to edge points should be updated to proper value according to
the program schemes. Take the nodes at the left edge of the array
($V_{i,1}, 1\leq i \leq n$) as an example. All of the coefficients
related to these nodes are located in $A1$, which should be set as

    \begin{equation}
    A1(k,k) = \left\{
    \begin{array}{ll}
    -(R_l^{-1}+R_{i,1}^{-1})   & \text{if } floating\\
    -(R_v^{-1}+R_l^{-1}+R_{i,1}^{-1})& \text{if } activated
    \end{array} \right.
    \end{equation}
    where $k=(n-1)i+1$ for $1\leq i\leq m$.

Similar procedures can be performed to initiate the coefficients
related to other edges.

\item The constant term $C$ is a $2mn{\times}1$ vector.
    Equation(\ref{equ:KCL1})-(\ref{equ:KCL8}) show that only KCL
    equations of the activated points have constant terms. Therefore,
    only the following elements in $C$ may have non-zero value:
    $C((i-1)n+1)$, $C(in)$, $C(mn+i)$ and $C((2m-1)n+i)$ for
    $i=1,2...m$, corresponding to the nodes at the four edges
    respectively. If the wordline is driven from one side, only nodes
    $V_{i,1}$ and $V_{1,j}$ can be activated. Therefore, the constant
    corresponding to these nodes can be defined as
\begin{equation}
    C((i-1)n+1) = \left\{
    \begin{array}{ll}
    0   & \text{if } floating\\
    -R_v^{-1}V_{Wi}& \text{if } activated
    \end{array} \right.
    \end{equation}
for wordline layer and
\begin{equation}
    C(mn+j) = \left\{
    \begin{array}{ll}
    0   & \text{if } floating\\
    -R_s^{-1}V_{Bj}& \text{if } activated
    \end{array} \right.
    \end{equation}
for bitline layer, with $1\leq i\leq m$, and $1\leq j\leq n$.
\end{enumerate}
